{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "naN9CKGafLWX"
      },
      "source": [
        "Copyright 2019 Google LLC\n",
        "\n",
        "Licensed under the Apache License, Version 2.0 (the \"License\");\n",
        "you may not use this file except in compliance with the License.\n",
        "You may obtain a copy of the License at\n",
        "\n",
        "    https://www.apache.org/licenses/LICENSE-2.0\n",
        "\n",
        "Unless required by applicable law or agreed to in writing, software\n",
        "distributed under the License is distributed on an \"AS IS\" BASIS,\n",
        "WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
        "See the License for the specific language governing permissions and\n",
        "limitations under the License."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "pHqqZAKjRjkn"
      },
      "outputs": [],
      "source": [
        "from __future__ import absolute_import\n",
        "from __future__ import division\n",
        "from __future__ import print_function\n",
        "\n",
        "import itertools\n",
        "import numpy\n",
        "import tensorflow as tf\n",
        "\n",
        "\n",
        "### Lie Algebra definitions for Spin(8), SU(8), E7.\n",
        "\n",
        "def permutation_sign(p):\n",
        "  q = [x for x in p]  # Copy to list.\n",
        "  parity = 1\n",
        "  for n in range(len(p)):\n",
        "    while n != q[n]:\n",
        "      qn = q[n]\n",
        "      q[n], q[qn] = q[qn], q[n]  # Swap to make q[qn] = qn.\n",
        "      parity = -parity\n",
        "  return parity\n",
        "\n",
        "\n",
        "class Spin8(object):\n",
        "  \"\"\"Container class for Spin(8) tensor invariants.\"\"\"\n",
        "\n",
        "  def __init__(self):\n",
        "    r8 = range(8)\n",
        "    self.gamma_vsc = gamma_vsc = self._get_gamma_vsc()\n",
        "    #\n",
        "    # The gamma^{ab}_{alpha beta} tensor that translates between antisymmetric\n",
        "    # matrices over vectors [ij] and antisymmetric matrices over spinors [sS].\n",
        "    # (These expressions could be shortened, but we want to keep\n",
        "    # symmetric-looking formulas.)\n",
        "    self.gamma_vvss = 0.5 * (\n",
        "        numpy.einsum('isc,jSc-\u003eijsS', gamma_vsc, gamma_vsc) -\n",
        "        numpy.einsum('jsc,iSc-\u003eijsS', gamma_vsc, gamma_vsc))\n",
        "    # The gamma^{ab}_{alpha* beta*} tensor that translates between antisymmetric\n",
        "    # matrices over vectors [ij] and antisymmetric matrices over cospinors [cC].\n",
        "    self.gamma_vvcc = 0.5 * (\n",
        "        numpy.einsum('isc,jsC-\u003eijcC', gamma_vsc, gamma_vsc) -\n",
        "        numpy.einsum('jsc,isC-\u003eijcC', gamma_vsc, gamma_vsc))\n",
        "    #\n",
        "    # The gamma^{ijkl}_{alpha beta} tensor that translates between antisymmetric\n",
        "    # 4-forms [ijkl] and symmetric traceless matrices over the spinors (sS).\n",
        "    g_ijsS = numpy.einsum('isc,jSc-\u003eijsS', self.gamma_vsc, self.gamma_vsc)\n",
        "    g_ijcC = numpy.einsum('isc,jsC-\u003eijcC', self.gamma_vsc, self.gamma_vsc)\n",
        "    g_ijklsS = numpy.einsum('ijst,kltS-\u003eijklsS', g_ijsS, g_ijsS)\n",
        "    g_ijklcC = numpy.einsum('ijcd,kldC-\u003eijklcC', g_ijcC, g_ijcC)\n",
        "    gamma_vvvvss = numpy.zeros([8] * 6)\n",
        "    gamma_vvvvcc = numpy.zeros([8] * 6)\n",
        "    for perm in itertools.permutations(range(4)):\n",
        "      perm_ijkl = ''.join('ijkl'[p] for p in perm)\n",
        "      sign = permutation_sign(perm)\n",
        "      gamma_vvvvss += sign * numpy.einsum(perm_ijkl + 'sS-\u003eijklsS', g_ijklsS)\n",
        "      gamma_vvvvcc += sign * numpy.einsum(perm_ijkl + 'cC-\u003eijklcC', g_ijklcC)\n",
        "    self.gamma_vvvvss = gamma_vvvvss / 24.0\n",
        "    self.gamma_vvvvcc = gamma_vvvvcc / 24.0\n",
        "\n",
        "  def _get_gamma_vsc(self):\n",
        "    \"\"\"Computes SO(8) gamma-matrices.\"\"\"\n",
        "    # Conventions match Green, Schwarz, Witten's, but with index-counting\n",
        "    # starting at zero.\n",
        "    entries = (\n",
        "        \"007+ 016- 025- 034+ 043- 052+ 061+ 070- \"\n",
        "        \"101+ 110- 123- 132+ 145+ 154- 167- 176+ \"\n",
        "        \"204+ 215- 226+ 237- 240- 251+ 262- 273+ \"\n",
        "        \"302+ 313+ 320- 331- 346- 357- 364+ 375+ \"\n",
        "        \"403+ 412- 421+ 430- 447+ 456- 465+ 474- \"\n",
        "        \"505+ 514+ 527+ 536+ 541- 550- 563- 572- \"\n",
        "        \"606+ 617+ 624- 635- 642+ 653+ 660- 671- \"\n",
        "        \"700+ 711+ 722+ 733+ 744+ 755+ 766+ 777+\")\n",
        "    ret = numpy.zeros([8, 8, 8])\n",
        "    for ijkc in entries.split():\n",
        "     ijk = tuple(map(int, ijkc[:-1]))\n",
        "     ret[ijk] = +1 if ijkc[-1] == '+' else -1\n",
        "    return ret\n",
        "\n",
        "\n",
        "class SU8(object):\n",
        "  \"\"\"Container class for su(8) tensor invariants.\"\"\"\n",
        "\n",
        "  def __init__(self):\n",
        "    # Tensor that translates between adjoint indices 'a' and\n",
        "    # (vector) x (vector) indices 'ij'\n",
        "    ij_map = [(i, j) for i in range(8) for j in range(8) if i \u003c j]\n",
        "    #\n",
        "    # We also need the mapping between 8 x 8 and 35 representations, using\n",
        "    # common conventions for a basis of the 35-representation, and likewise\n",
        "    # for 8 x 8 and 28.\n",
        "    m_35_8_8 = numpy.zeros([35, 8, 8], dtype=numpy.complex128)\n",
        "    m_28_8_8 = numpy.zeros([28, 8, 8], dtype=numpy.complex128)\n",
        "    for n in range(7):\n",
        "      m_35_8_8[n, n, n] = +1.0\n",
        "      m_35_8_8[n, n + 1, n + 1] = -1.0\n",
        "    for a, (m, n) in enumerate(ij_map):\n",
        "      m_35_8_8[a + 7, m, n] = m_35_8_8[a + 7, n, m] = 1.0\n",
        "      m_28_8_8[a, m, n] = 1.0\n",
        "      m_28_8_8[a, n, m] = -1.0\n",
        "    #\n",
        "    # The su8 'Generator Matrices'.\n",
        "    t_aij = numpy.zeros([63, 8, 8], dtype=numpy.complex128)\n",
        "    t_aij[:35, :, :] = 1.0j * m_35_8_8\n",
        "    for a, (i, j) in enumerate(ij_map):\n",
        "      t_aij[a + 35, i, j] = -1.0\n",
        "      t_aij[a + 35, j, i] = 1.0\n",
        "    self.ij_map = ij_map\n",
        "    self.m_35_8_8 = m_35_8_8\n",
        "    self.m_28_8_8 = m_28_8_8\n",
        "    self.t_aij = t_aij\n",
        "\n",
        "\n",
        "class E7(object):\n",
        "  \"\"\"Container class for e7 tensor invariants.\"\"\"\n",
        "\n",
        "  def __init__(self, spin8, su8):\n",
        "    self._spin8 = spin8\n",
        "    self._su8 = su8\n",
        "    ij_map = su8.ij_map\n",
        "    t_a_ij_kl = numpy.zeros([133, 56, 56], dtype=numpy.complex128)\n",
        "    for a in range(35):\n",
        "      # In order to speed up initialization, we are splitting X x Y x Z\n",
        "      # tensor products, which numpy.einsum() would compute without generating\n",
        "      # intermediate products, into sequences of X x Y tensor products.\n",
        "      t_a_ij_kl[:35, 28:, :28] = (1 / 8.0) * (\n",
        "          numpy.einsum('qIkl,Kkl-\u003eqIK',\n",
        "                       numpy.einsum(\n",
        "                           'ijklq,Iij-\u003eqIkl',\n",
        "                           numpy.einsum('ijklsS,qsS-\u003eijklq',\n",
        "                                        spin8.gamma_vvvvss, su8.m_35_8_8),\n",
        "                           su8.m_28_8_8),\n",
        "                       su8.m_28_8_8))\n",
        "      t_a_ij_kl[:35, :28, 28:] = (1 / 8.0) * (\n",
        "          numpy.einsum('qIkl,Kkl-\u003eqIK',\n",
        "                       numpy.einsum(\n",
        "                           'ijklq,Iij-\u003eqIkl',\n",
        "                           numpy.einsum('ijklsS,qsS-\u003eijklq',\n",
        "                                        spin8.gamma_vvvvss, su8.m_35_8_8),\n",
        "                           su8.m_28_8_8),\n",
        "                       su8.m_28_8_8))\n",
        "      #\n",
        "      t_a_ij_kl[35:70, 28:, :28] = (1.0j / 8.0) * (\n",
        "          numpy.einsum('qIkl,Kkl-\u003eqIK',\n",
        "                       numpy.einsum(\n",
        "                           'ijklq,Iij-\u003eqIkl',\n",
        "                           numpy.einsum('ijklcC,qcC-\u003eijklq',\n",
        "                                        spin8.gamma_vvvvcc, su8.m_35_8_8),\n",
        "                           su8.m_28_8_8),\n",
        "                       su8.m_28_8_8))\n",
        "      t_a_ij_kl[35:70, :28, 28:] = (-1.0j / 8.0) * (\n",
        "          numpy.einsum('qIkl,Kkl-\u003eqIK',\n",
        "                       numpy.einsum(\n",
        "                           'ijklq,Iij-\u003eqIkl',\n",
        "                           numpy.einsum('ijklcC,qcC-\u003eijklq',\n",
        "                                        spin8.gamma_vvvvcc, su8.m_35_8_8),\n",
        "                           su8.m_28_8_8),\n",
        "                       su8.m_28_8_8))\n",
        "      #\n",
        "      # We need to find the action of the su(8) algebra on the\n",
        "      # 28-representation.\n",
        "      su8_28 = 2 * (\n",
        "          numpy.einsum(\n",
        "              'aIjn,Jjn-\u003eaIJ',\n",
        "              numpy.einsum(\n",
        "                  'aimjn,Iim-\u003eaIjn',\n",
        "                  numpy.einsum(\n",
        "                      'aij,mn-\u003eaimjn',\n",
        "                      su8.t_aij,\n",
        "                      numpy.eye(8, dtype=numpy.complex128)),\n",
        "                  su8.m_28_8_8),\n",
        "              su8.m_28_8_8))\n",
        "      t_a_ij_kl[70:, :28, :28] = su8_28\n",
        "      t_a_ij_kl[70:, 28:, 28:] = su8_28.conjugate()\n",
        "      self.t_a_ij_kl = t_a_ij_kl\n",
        "\n",
        "\n",
        "spin8 = Spin8()\n",
        "su8 = SU8()\n",
        "e7 = E7(spin8, su8)\n",
        "\n",
        "def get_proj_35_8888(want_selfdual=True):\n",
        "  \"\"\"Computes the (35, 8, 8, 8, 8)-projector to the (anti)self-dual 4-forms.\"\"\"\n",
        "  # We first need some basis for the 35 self-dual 4-forms.\n",
        "  # Our convention is that we lexicographically list those 8-choose-4\n",
        "  # combinations that contain the index 0.\n",
        "  sign_selfdual = 1 if want_selfdual else -1\n",
        "  ret = numpy.zeros([35, 8, 8, 8, 8], dtype=numpy.float64)\n",
        "  #\n",
        "  def get_selfdual(ijkl):\n",
        "    mnpq = tuple(n for n in range(8) if n not in ijkl)\n",
        "    return (sign_selfdual * permutation_sign(ijkl + mnpq),\n",
        "            ijkl, mnpq)\n",
        "  selfduals = [get_selfdual(ijkl)\n",
        "               for ijkl in itertools.combinations(range(8), 4)\n",
        "               if 0 in ijkl]\n",
        "  for num_sd, (sign_sd, ijkl, mnpq) in enumerate(selfduals):\n",
        "    for abcd in itertools.permutations(range(4)):\n",
        "      sign_abcd = permutation_sign(abcd)\n",
        "      ret[num_sd,\n",
        "          ijkl[abcd[0]],\n",
        "          ijkl[abcd[1]],\n",
        "          ijkl[abcd[2]],\n",
        "          ijkl[abcd[3]]] = sign_abcd\n",
        "      ret[num_sd,\n",
        "          mnpq[abcd[0]],\n",
        "          mnpq[abcd[1]],\n",
        "          mnpq[abcd[2]],\n",
        "          mnpq[abcd[3]]] = sign_abcd * sign_sd\n",
        "  return ret / 24.0\n",
        "\n",
        "\n",
        "### Supergravity.\n",
        "\n",
        "def tf_so8_sugra_stationarity(t_a1, t_a2):\n",
        "  \"\"\"Computes the stationarity-condition tensor.\"\"\"\n",
        "  # See: https://arxiv.org/pdf/1302.6219.pdf, text after (3.2).\n",
        "  t_x0 = (\n",
        "      +4.0 * tf.einsum('mi,mjkl-\u003eijkl', t_a1, t_a2)\n",
        "      -3.0 * tf.einsum('mnij,nklm-\u003eijkl', t_a2, t_a2))\n",
        "  t_x0_real = tf.real(t_x0)\n",
        "  t_x0_imag = tf.imag(t_x0)\n",
        "  tc_sd = tf.constant(get_proj_35_8888(True))\n",
        "  tc_asd = tf.constant(get_proj_35_8888(False))\n",
        "  t_x_real_sd = tf.einsum('aijkl,ijkl-\u003ea', tc_sd, t_x0_real)\n",
        "  t_x_imag_asd = tf.einsum('aijkl,ijkl-\u003ea', tc_asd, t_x0_imag)\n",
        "  return (tf.einsum('a,a-\u003e', t_x_real_sd, t_x_real_sd) +\n",
        "          tf.einsum('a,a-\u003e', t_x_imag_asd, t_x_imag_asd))\n",
        "\n",
        "\n",
        "def tf_so8_sugra_potential(t_v70):\n",
        "  \"\"\"Returns dict with key tensors from the SUGRA potential's TF graph.\"\"\"\n",
        "  tc_28_8_8 = tf.constant(su8.m_28_8_8)\n",
        "  t_e7_generator_v70 = tf.einsum(\n",
        "      'v,vIJ-\u003eJI',\n",
        "      tf.complex(t_v70, tf.constant([0.0] * 70, dtype=tf.float64)),\n",
        "      tf.constant(e7.t_a_ij_kl[:70, :, :], dtype=tf.complex128))\n",
        "  t_complex_vielbein = tf.linalg.expm(t_e7_generator_v70)\n",
        "  def expand_ijkl(t_ab):\n",
        "    return 0.5 * tf.einsum(\n",
        "        'ijB,BIJ-\u003eijIJ',\n",
        "        tf.einsum('AB,Aij-\u003eijB', t_ab, tc_28_8_8),\n",
        "        tc_28_8_8)\n",
        "  #\n",
        "  t_u_ijIJ = expand_ijkl(t_complex_vielbein[:28, :28])\n",
        "  t_u_klKL = expand_ijkl(t_complex_vielbein[28:, 28:])\n",
        "  t_v_ijKL = expand_ijkl(t_complex_vielbein[:28, 28:])\n",
        "  t_v_klIJ = expand_ijkl(t_complex_vielbein[28:, :28])\n",
        "  #\n",
        "  t_uv = t_u_klKL + t_v_klIJ\n",
        "  t_uuvv = (tf.einsum('lmJK,kmKI-\u003elkIJ', t_u_ijIJ, t_u_klKL) -\n",
        "            tf.einsum('lmJK,kmKI-\u003elkIJ', t_v_ijKL, t_v_klIJ))\n",
        "  t_T = tf.einsum('ijIJ,lkIJ-\u003elkij', t_uv, t_uuvv)\n",
        "  t_A1 = (-4.0 / 21.0) * tf.trace(tf.einsum('mijn-\u003eijmn', t_T))\n",
        "  t_A2 = (-4.0 / (3 * 3)) * (\n",
        "      # Antisymmetrize in last 3 indices, taking into account antisymmetry\n",
        "      # in last two indices.\n",
        "      t_T\n",
        "      + tf.einsum('lijk-\u003eljki', t_T)\n",
        "      + tf.einsum('lijk-\u003elkij', t_T))\n",
        "  t_A1_real = tf.real(t_A1)\n",
        "  t_A1_imag = tf.imag(t_A1)\n",
        "  t_A2_real = tf.real(t_A2)\n",
        "  t_A2_imag = tf.imag(t_A2)\n",
        "  t_A1_potential = (-3.0 / 4) * (\n",
        "      tf.einsum('ij,ij-\u003e', t_A1_real, t_A1_real) +\n",
        "      tf.einsum('ij,ij-\u003e', t_A1_imag, t_A1_imag))\n",
        "  t_A2_potential = (1.0 / 24) * (\n",
        "      tf.einsum('ijkl,ijkl-\u003e', t_A2_real, t_A2_real) +\n",
        "      tf.einsum('ijkl,ijkl-\u003e', t_A2_imag, t_A2_imag))\n",
        "  t_potential = t_A1_potential + t_A2_potential\n",
        "  #\n",
        "  return dict(v70=t_v70,\n",
        "              vielbein=t_complex_vielbein,\n",
        "              tee_tensor=t_T,\n",
        "              a1=t_A1,\n",
        "              a2=t_A2,\n",
        "              potential=t_potential)\n",
        "\n",
        "\n",
        "def call_with_critical_point_scanner(f, *args):\n",
        "  \"\"\"Calls f(scanner, *args) in TensorFlow session-context.\n",
        "\n",
        "  Here, `scanner` will be a function with signature\n",
        "  scanner(seed:int, scale:float) -\u003e (potential, stationarity, pos_vector).\n",
        "\n",
        "  The function `scanner` can only perform a scan when called from within\n",
        "  the TF session-context that is set up by this function.\n",
        "  \"\"\"\n",
        "  graph = tf.Graph()\n",
        "  with graph.as_default():\n",
        "    t_input = tf.placeholder(tf.float64, shape=[70])\n",
        "    t_v70 = tf.Variable(initial_value=numpy.zeros([70]),\n",
        "                        trainable=True,\n",
        "                        dtype=tf.float64)\n",
        "    op_assign_input = tf.assign(t_v70, t_input)\n",
        "    d = tf_so8_sugra_potential(t_v70)\n",
        "    t_potential = d['potential']\n",
        "    t_stationarity = tf_so8_sugra_stationarity(d['a1'], d['a2'])\n",
        "    opt = tf.contrib.opt.ScipyOptimizerInterface(tf.asinh(t_stationarity),\n",
        "                                                 options=dict(maxiter=500))\n",
        "    with tf.Session() as sess:\n",
        "      sess.run([tf.global_variables_initializer()])\n",
        "      def scanner(seed, scale):\n",
        "        rng = numpy.random.RandomState(seed)\n",
        "        v70 = rng.normal(scale=scale, size=[70])\n",
        "        sess.run([op_assign_input], feed_dict={t_input: v70})\n",
        "        opt.minimize(session=sess)\n",
        "        n_ret = sess.run([t_potential, t_stationarity, t_v70])\n",
        "        return n_ret\n",
        "      return f(scanner, *args)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "colab": {
          "height": 364
        },
        "colab_type": "code",
        "executionInfo": {
          "elapsed": 61136,
          "status": "ok",
          "timestamp": 1559567551803,
          "user": {
            "displayName": "",
            "photoUrl": "",
            "userId": ""
          },
          "user_tz": -120
        },
        "id": "paC8SebISgVX",
        "outputId": "869e5c4f-eb6a-44db-875c-52993796a66a"
      },
      "outputs": [],
      "source": [
        "### Demo. Scans for 20 critical points.\n",
        "\n",
        "\n",
        "def scan_many(scanner, num_to_produce):\n",
        "  ret = []\n",
        "  seed = 0\n",
        "  while len(ret) \u003c num_to_produce:\n",
        "    seed += 1\n",
        "    scanned = scanner(seed, scale=2.0)\n",
        "    potential, stationarity, v70 = scanned\n",
        "    if stationarity \u003c 1e-3:\n",
        "      ret.append(scanned)\n",
        "      print('%3d / %3d: V=%.6f stationarity=%.3g' % (\n",
        "          len(ret), num_to_produce, potential, stationarity))\n",
        "  return ret\n",
        "\n",
        "\n",
        "scanned = call_with_critical_point_scanner(scan_many, 20)\n"
      ]
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "name": "so8_supergravity.ipynb",
      "provenance": [],
      "version": "0.3.2"
    },
    "kernelspec": {
      "display_name": "Python 2",
      "name": "python2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
